SAGA-HE-117-97 (|hep-ph/97Uy22U|) 



October 28, 1997 



Numerical solution of evolution equations 
for polarized structure functions 

M. Hirai, S. Kumano, and M. Miyama * 
Department of Physics, Saga University, Saga 840, Japan 

ABSTRACT 

We investigate numerical solution of Dokshitzer-Gribov-Lipatov-Altarelli-Parisi 
(DGLAP) evolution equations for longitudinally polarized structure functions. Fla- 
vor nonsinglet and singlet equations with next-to-leading-order as corrections are stud- 
ied. A brute-force method is employed. Dividing the variables x and into small 
steps, we simply solve the integrodifferential equations. Numerical results indicate that 
accuracy is better than 1% in the region 10~^ < x < 0.8 if more than two-hundred 

steps and more than one-thousand x steps are taken. Our evolution results are 
compared with polarized experimental data of the spin asjTiimetry Ai by the SLAC- 
E130, SLAC-E143, EMC, and SMC collaborations. The comparison indicates that we 
cannot assume Ai is independent of Q^. We provide a FORTRAN program for the 
evolution and devolution of polarized nonsinglet-quark, singlet-quark, Ag, -t- Agj, and 
gluon distributions (and corresponding structure functions). 
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Program Summary 



Title of program: BFPl 

Computer: AlphaServer 2100 4/200; Installation: The Research Center for Nuclear 
Physics in Osaka 

Operating system: Open VMS V6.1 
Programming language used: FORTRAN 77 
Peripherals used: Laser printer 

No. of lines in distributed program, including test data, etc.: 1617 

Keywords: Structure function, polarized parton distribution, evolution, numerical 
solution. 

Nature of physical problem 

This program solves DGLAP evolution equations with or without next-to-leading- 
order as effects for longitudinally polarized parton distributions. The evolved distri- 
butions could be convoluted with coefficient functions for calculating the structure 
function gi. Both flavor-nonsinglet and singlet cases are provided, so that the distri- 
butions, xAq^^g, xAq^, xAq^ = xAqi + xAqi (i=quark flavor), xAg, xgi^^, xgig, and 
xg^^ can be obtained. 

Method of solution 

The DGLAP integrodifferential equations are simply solved by a brute-force method. 
We divide the variables x and into very small steps, and the integrodifferential 
equations are solved step by step. 

Restrictions of the program 

This program is used for calculating evolution of a longitudinally polarized flavor- 
nonsinglet-quark, singlet-quark, Aq^, and gluon distributions in the leading order or 
in the next-to-leading order of a^. evolution equations are the DGLAP equations. 
The double precision arithmetic is used. The renormalization scheme is the modified 
minimal subtraction scheme (MS). A user provides the initial distribution as a sub- 
routine or as a data file. Examples are explained in section ^. Then, the user inputs 
fifteen parameters in section 

Typical running time 

Approximately seven minutes on AlphaServer 2100 4/200 in the nonsinglet case, sixty 
minutes in the singlet-quark evolution. 
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LONG WRITE-UP 



1 Introduction 

Spin structure of the proton has been investigated by polarized lepton-proton scatter- 
ing. In particular, the structure function gi provides us information on the probability 
to find a quark polarized along the direction of longitudinally polarized proton spin 
minus the one polarized oppositely. Experimentally, it is measured by a spin asym- 
metry in the lepton-proton scattering. The gi is a function of two kinematical vari- 
ables and x = Q'^/(2p ■ q), where p is the proton momentum and is given by 
the four- momentum transfer g as = — Q'^- The x dependence is associated with 
nonperturbative physics, so that x distributions of polarized parton distributions are 
obtained by fitting various experimental data. At this stage, the only way to study 
the parton x distributions theoretically is to use phenomenological quark models. On 
the other hand, the dependence is well predicted in perturbative QCD. In general, 
structure functions are almost independent of the scale Q^, which is called Bjorken 
scaling. However, even though it is small, scaling-violation phenomena are observed. 
The experimental data support perturbative QCD predictions. 

As a way of describing the scaling violation, Dokshitzer-Gribov-Lipatov-Altarelli- 
Parisi (DGLAP) equations ^ are usually used. They are coupled integrodifferential 
equations. Because the DGLAP equations are frequently used in theoretical and ex- 
perimental studies, it is worth while creating a program to solve them accurately. We 
have been studying this topic last several years in a Laguerre-polynomial method [Q] 
and in a brute-force method p, If. Here, we extend the studies to the spin- dependent 
case. Fortunately, the next-to-leading-order (NLO) splitting functions are calculated 
recently for polarized parton distributions 0. It enables us to investigate the details 
of NLO effects on the evolution. We create a FORTRAN program for obtaining 
numerical solution of the polarized DGLAP equations with the NLO effects. With 
our program, evolution of nonsinglet quark, singlet quark, and gluon distributions 
can be obtained. Furthermore, NLO coefficient functions could be convoluted with the 
evolved distributions in our program for calculating the gi structure function. It is 
very useful for investigating polarized parton distributions. In fact, our program in a 
slightly modified form is used for investigating optimum polarized parton distributions 
by fitting experimental data 

We explain the details of our studies in the following. In section ^, the polarized 
DGLAP evolution equations are explained. Then, our numerical solution method is 
described in section ^. Input parameters and input distributions are discussed in section 
U, and essential subroutines in the program are explained in section |^. Numerical results 
and their comparison with experimental data are discussed in section ^. Summary is 
given in section |^. Explicit forms of splitting functions and the coefficient functions 
are listed in Appendices. 
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2 evolution equations 

We use the DGLAP equations ||l|] for studying the evolution. Both the leading 
order (LO) and the next-to-leading order (NLO) cases can be handled by the DGLAP 
equations. NLO effects are included in the running coupling constant ^^(Q^) and in 
the splitting functions APjj (x). Here, the evolution of polarized parton distributions 
Ag = ^-f — and A^f = (7+1 — is investigated. 
First, the nonsinglet DGLAP equation is given by 



2^ 



Ag^,(x,Q^) = AP,±,;v5(x) ® Ag^,(x,Q^) , (2.1) 

where Ag^g(x, Q^) is a polarized nonsinglet parton distribution, APg± tvs is the polar- 
ized nonsinglet splitting function, and the convolution is defined by 



X 



/(a;)®^(x)= / -J^fl-\g{y) . (2.2) 

The notation in the splitting function indicates a Ag"*" = Ag + Ag or Ag~ = Ag — Ag 
distribution type, which is explained in Appendix A. Instead of Q^, it is more convenient 
to use the variable t defined by 



2 , 



«s(g^) 



(2.3) 



where /9o is defined in Appendix A. The parton distribution and the splitting function 
multiplied by x 

fix) = xfix) (2.4) 

satisfy the same integrodifferential equation. Therefore, we rewrite the evolution equa- 
tion as 

d ^ 

— Ag^^ (x, t) = APg±^Nsix) ® Ag^^ (x, t) . (2.5) 

Next, the singlet evolution is more complicated than the nonsinglet one due to gluon 
participation in the evolution. The singlet quark distribution is defined by Aqs{x,t) = 
X Aqf where i is the flavor, and Aljix, t) is the gluon distribution. The singlet 
case is given by 

d_ f Ag,(x,t) \ ^ ( APg,{x,t) APg,ix,t) \ f Ag;(x,t) \ 

dt \ Ag{x,t) J \ AP,,{x,t) APgg{x,t) i V ^9{x,t) ) ' ^ ' ' 
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If each flavor evolution is necessary, another equation has to be used in addition to Eq. 



. Using the properties of the splitting functions, AP +^ + = 6ij AP^+ +2Ci?TRFc 



q^,Ns 
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5(1) 
d_ 



and ^P^j+g = ^Pq+g (independent of j), we have the evolution equation 



(2.7) 

In this case, three coupled integrodifferential equations in Eqs. ( P^ ) and ( |2.7| ) should 
be solved simultaneously. All the necessary splitting functions AP and F are listed in 
Appendix A. 

We discuss the NLO effects in the evolution equations. The NLO contributions are 
included in the running coupling constant as{t), in the splitting functions APij{z), and 
in the coefficient functions. Once the NLO corrections are included in the evolution, 
the renormalization scheme has to be specified. We use the MS scheme throughout 
this paper. 

The running coupling constant in the leading order (LO) is given by 



and the one in the next-to-leading order (NLO) is 



a 



NLO 



Air 



/?oln(gVA2 



1 - 



Aln{ln(gVA^)} 
/?o'ln(gVA2) 



(2. 



(2.9) 



The constants /3o and (3i are defined in Appendix A. The splitting functions have the 
perturbative-expansion form 



AP,,(x) 



A/f)(x) + ^Aif)(x) 



(2.10) 



The second term is the NLO contributions to the splitting functions. The functions 
are expressed by the ones multiplied by x (AP). Changing the variable to t in the 
DGLAP equations, we have the splitting functions 



AP,,(x) 



where the function APjj(x) is 



AR^j{x) 



AP?^(x) + 



27r 



ARij{x) , 



^ApWfx) 
2/?o ^ ' 



(2.11) 



f2.12) 



The second term in Eq. ( |2.12|) appears because of the transformation from to t. To 
be precise, the splitting functions APj^ in Eq. ( |2.11| ) should be denoted, for example. 



AP^- because it is different from xAP.j = xAP^f + {as/2n)xAPl;f which is Eq. (|2A0| ) 

multiphed by x. (Note the definition of / is / = xf .) However, we use the expression 
without the prime throughout this paper for simphcity. We have the same expression 
APns = APj^l + as/{2TT)ARNs in the nonsinfflet case. 

Next, we discuss the spin-dependent structure function gi. We have discussed how 
the evolution of quark and gluon distributions is described. However, these parton 
distributions are not observed directly in experiments. The gi could be measured in 
the polarized lepton-proton scattering with information on the unpolarized structure 
function Fi. In the LO case, it is given in parton model as 

Nf 

^^?i(^) = ^E^'^^''(^)- (2.13) 

i 

The NLO effects in the structure function are included in the coefficient functions and 
also in the quark and gluon distributions. In the NLO case, the quark distributions 
should be convoluted with a coefficient function. Furthermore, an additional gluon 
correction term should be taken into account: 

xgi{x) = -J2e^, \^AC,ix)^Aq + {x)+ACgix)^Ag{x)j , (2.14) 

i 

where ACq/x and ACg/x are the quark and gluon coefficient functions in Appendix B. 
The calculation procedure for the gi evolution is the following. First, the evolution 
of the quark and gluon distributions is calculated. Then, the structure function gi at 
is evaluated by using Eq. ( p.l3| ) or ( p.l4|) . An example is explained in section |6l2 . 
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3 Brute-force method 



Among various methods of solving the DGLAP equations, we decide to employ a brute- 
force method. So far, we have investigated two methods, the Laguerre-polynomial [0] 
and the brute-force § methods, in the spin-independent case. The Laguerre method 
has an advantage of computing time. However, numerical accuracy becomes slightly 
worse in the nonsinglet case at small x. As far as we studied, the situation is better in 
the polarized distributions. However, we find a tendency that the accuracy is slightly 
worse at small x. In light of future HERA spin physics, our program should be accu- 
rately enough even at x = 10^^. Therefore, we consider that the brute-force method is 
safer for getting accurate results in the wide x range. The variables x and t are divided 
into small steps, then integration and differentiation are defined by 

dfjx) _ f{Xm+l) - fjXm) , ^ 



f{x) dx=^ SXmfiXm) ■ (3.2) 
m=l 



With these replacements, the evolution equations could be solved rather easily. For 
example, Eq. (p.5|) is written in the following form: 



^QNsi^k,tj+i) = Aq^g{xk,tj) + Stj -^APg±^Ns ( — I Aq^g{xm,tj) . (3.3) 

, Xm \Xm I 

m=k ' 

First, the evolution from = to ti = 5t is calculated in the above equation by 
providing the initial distribution Ag^g(xm, ij=i)- Repeating this step A^t — 1 times, we 
obtain the final distribution at tjq^. Accurate results cannot be obtained at small x if 
the linear x step is taken [6x = l/N^). Therefore, the logarithmic-x step 6{logiQx) = 
\logioXmin\/Nx is taken in our evolution calculations. 

The same method can be applied to the singlet and Aq^' evolution equations in 
Eqs. (p^ ) and (|2.7| ). These equations are written in the brute- force method as: 

Ag. + (xfc,tj+i) =Aq^^{xk,tj) + 6tj 22^^^Pg+,NS I — I ^Q^^{xm,tj) 

m—k "^"^ \Xm / 

+ 6t, Y.^-^2CpTnF,, Ag,(x™,t,) 

m=k 

+ 6t, ^ ^ AP^.^ (^] Ag{x^, t,) , (3.4a) 

m—k \Xm / 
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5Xr 



m=k 



Aq^{Xk,tj+i) ^Aq^{Xk,tj) + 5tj ^^APqq ( — ) Aqg{Xm,tj) 



Sx. 



m=k 



+ 5]^AP,,(^) Ag(x^,t,), 



(3.4b) 



5Xr. 



m=k 



Ag{xk,tj+i) =Ag{xk,tj) + Stj APgq f — ) Aqg{xm,t_ 



+ 5^ ^ AP,, ) Ag{x^,t,). 



m=k 



(3.4c) 



These are coupled equations. However, they can be solved by providing the initial 
distributions, Aq^{xjn-itj=Q) , Aqg{xm,tj=o), and Ag{xm,tj=o), and by repeating the 
evolution step Nt times. 

We should be careful in handling 1/(1 — a;)+ terms in the splitting functions. They 
are given in our method by 



dx' 



(1 - x'U 



f{Xm) - /(I) 



m=k 



Xv 



+ /(l)ln(l-a;fe) . 



In the same way, the integral with [ln(l — x)/{\ — a;)]+ becomes 



dx'f{x') 



ln(l - x') 
l-x' 



E 6x^ [f{xm) - /( 



+ m=k 
1 



I- Xr, 



+ 2/(1) 



(3.5) 



(3.6) 
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4 Description of input parameters and input 
distribution 



For running the FORTRAN-77 program BFPl, a user should supply fifteen input 
parameters from the file #10. In addition, an input distribution (s) should be given 
in a function subroutine(s) in the end of the FORTRAN program or in an input data 
file(s), #13, #14, and/or #15. The initial distribution (s) could be written in the 
output file #12. Evolution results are written in the output file #11. We explain the 
input parameters and the input distributions in the following. 

4.1 Input parameters 



There are fifteen parameters. Numerical values of the parameters should be supplied 
in the file #10, then these arc read in the main program. 



lOUT 


1 


write x and xAq^g{x) [or xgi^,,{x)] at fixed (=Q2) in the file #11 


2 


write Q^ and xAq^g{Q'^) [or xgij^g{Q'"^)] at fixed x (=XX) in the file #11 


3 


X, xAgg(x) [or xgig{x)], and xAg{x) 


4 


Q^ xAq^iQ'^) [or xgig(Q^)], and xAg{Q'^) 


5 


X, xAq^{x) [or xg^-{x)], xAq^^{x) [or xgig{x)], and xAg{x) [or xgi^g{x)] 


6 


Q\ xAqtiQ-') [or (Q^)], xAq,{Q') [or xgi.iQ')], 
and xAg{Q'^) [or xgi^g{x)\ 


IREAD 


1 


give initial distribution (s) in function subroutine(s) 


2 


read initial distribution(s) from data file(s) 


INDIST 


1 


do not write initial distribution(s) 


2 


write initial distribution(s) in the file #12 


3 


write initial distribution(s) in the file #12 without calculating evolution 


lORDER 


1 


leading order (LO) in ag 


2 


u('xl-l()-l('a,dius order (NLO) 


ITYPE 


1 


structure function xgi{x,Q^) 


2 


quark distribution xAq{x, Q^) and gluon distribution xAg{x, Q^) 


IMORP 


1 


xAqi — xAqi type distribution 


2 


xAqi + xAqi type distribution 


Q02 


iniliMl (= Qq in GcV" ) ;il which an initial ditilribnlion iti snpplicd 




Q2 


Q'^ to which the distribution is evolved {Q'^ ^ Qq, Q'^ > Qq or Q'^ < Qq) 


DLAM 


QCD scale parameter Aqqd in GeV 


NF 


number of quark flavors 


XX 


X at which dependent distributions are written (I0UT=2, 4, or 6 case) 


NX 


number of x steps (NX<5000) 


NT 


number of t steps (NT<5000) 


NSTEP 


number of x steps or t steps for writing output distribution(s) 


NXMIN 


to(7io (minimum of x) [0 < min{x) = iq"^"^'" <^ XX] 
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The meaning of IREAD is explained in section 4.2. The structure function xQi 
is obtained by the convolution of the distributions, (l/2)[(4/9)x(A-u + Au + Ac + 
Ac) + (l/9)x(A(i + Ad + As + As)] and xAg in the four flavor case, with the corre- 
sponding coefficient functions. Practically, we use the expression {l/2)Y2-efxAqi = 
(l/2)a;[4Ag, - 3(Arf+ + As+)]/9, where Arf+ = Ad + Ad and As+ = As + As, instead 
of the above one in calculating xgi. It is explained in more detail in section |6.2| . The 
expression could be used only in the four and three flavor cases. In the five or six fiavor 
evolution, it should be slightly modified. The parameter IMORP indicates a plus or 
minus type distribution aix{Aqi±Aqi), where Oj are some constants. IM0RP=1 or 
2 should be taken in the nonsinglet evolution. In the singlet or Ag^^ case, IM0RP=2 
should be chosen. 

For example, if one would like to evolve an initial singlet-quark distribution xAq^ 
at Q2=4 GeV^ to the distribution at Q'^=200 GeV^ by the NLO DGLAP equations 
with Nf=4: and A=0.231 GeV, the input parameters could be I0UT=3, IREAD=1, 
INDIST=1, I0RDER=2, ITYPE=2, IM0RP=2, Q02=4.0, Q2=200.0, DLAM=0.231, 
NF=4, XX=0.0, NX=1000, NT=200, NSTEP=100, and NXMIN=-4. In this case, 
the input file #10 is the following: 

3, 1, 1, 2, 2, 2 

4.0, 200.0, 0.231, 4, 0.0, 1000, 200, 100, -4 . 

4.2 Input distributions supplied by function subroutines 
(IREAD=1) 

If IREAD=1 is chosen, an input distribution (s) at Ql should be supplied in the end of 
the FORTRAN program BFPl as a function subroutine (s). 

1) Nonsinglet case 

An initial nonsinglet- quark distribution at Qq should be given in QNSO(X) as a 
double precision function. As an example, the GS (set A) valence quark distribution 
xAu^ + xAd^ at Ql=4: GeV^ is given in the program BFPl. 

2) Singlet case 

An initial singlet-quark distribution at Ql should be supplied in QSO(X), and an 
initial gluon distribution in the nucleon should be in GO(X). The GS xAq^ = xAuy + 
xAdy + 6xAS and xAg distributions are given in the program. 

3) Aql' distribution case 

In addition to the above xAq^ and xAg distributions, the initial xAq^ distribution 
(and another fiavor distribution xAqj~) should be supplied. In calculating two-fiavor 
distributions simultaneously for obtaining xgi, two distributions (e.g. xAd~^ and a;As"'") 
should be supplied in the functions QIO(l,x) and QI0(2,x). If one needs only one-fiavor 
evolution, one may set QI0(2,x)=0. The GS distributions xAd'^ and xAs~^ are given 
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in our program as an example. 



4.3 Input distributions supplied by data files (IREAD=2) 

If IREAD=2 is chosen, an input distribution(s) at Ql should be supphed in a separate 
data file(s). 

1) Nonsinglet case (data file #13) 

An initial nonsinglet-quark distribution at Qq should be given in the data file #13 
as shown in the following example. 



0.000100 

0.000110 
0.000120 
0.000132 
0.000145 



0.009761 

0.010184 
0.010624 
0.011081 
0.011555 



1.000000 



0.000000 



The first column is the x values and the second one is the corresponding x/!^Uy + 

and at x—l.Q must be supplied. 



x/S.dy values. The data at x < Xmin 



2) Singlet case in the nucleon (data file #14) 

An initial singlet-quark distribution and a gluon distribution should be given in the 
data file #14 as shown in the following. 



0.000100 
0.000110 

0.000120 
0.000132 
0.000145 



0.001721 
0.001612 

0.001485 
0.001339 
0.001172 



0.006729 
0.007193 

0.007688 
0.008218 
0.008784 



1.000000 



0.000000 



0.000000 



The first column is the x values, the second is the xAq^ distribution, and the third 
is the x/S.g distribution. The data at x < x^m and at x—l.Q must be supplied. 

3) /S.ql distribution case (data file #15) 

An initial ^qf distribution, a singlet-quark distribution, and a gluon distribution 
should be given in the data file #15 as shown in the following. 
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0.000100 
0.000110 
0.000120 
0.000132 
0.000145 



-0.003637 
-0.003884 
-0.004147 
-0.004428 
-0.004727 



-0.002681 
-0.002859 
-0.003048 
-0.003249 
-0.003462 



0.001721 
0.001612 
0.001485 
0.001339 
0.001172 



0.006729 
0.007193 
0.007688 
0.008218 
0.008784 



1.000000 



0.000000 



0.000000 



0.000000 



0.000000 



The first column is the x values, the second is the x/S.q^ distribution, the third is 
another flavor distribution xAqj~, the fourth is xAq^, and the fifth is xAg. The data 
at a; < Xmin and at x=1.0 must be supplied. If one needs only one-flavor evolution, the 
second or third column values are set to 0. 



11 



5 Description of the program BFPl 



The major part of the program BFPl is essentially the same as the unpolarized version 
0. We do not explain each subroutine in this section. They are already well described 
in the previous publication, so that the interested reader may read Ref. for the 
detailed description of all the subroutines. 

The main program reads fifteen input parameters from the input file #10. Actual 
evolution calculations are done by calling the subroutine GETQNS in the nonsinglet 
case and GETQS in the singlet (or sAgj) case. Each evolution step is calculated by 
the subroutine QNSXT or GETQGX. Minor modifications are made in each flavor 
evolution (GETQS) so that two quark distributions [x!\ql and a;Ag^) are evolved 
simultaneously. A typical function of the GS parton distributions is given as the 
function GS(X,A,B,C,D,E,F). We explain the subroutines associated with the input 
distributions in the following. 



5.1 Function GS(X,A,B,C,D,E,F) 

This function calculates the GS parton distributions given by the parameters A, B, C, 
D, E, and F as ABx^{\ - x)^{\ + Ex + F^). 

5.2 Functions QNSO(X), QSO(X), GO(X), and QIO(I,X) 



The functions QNSO, QSO, and GO calculate an initial nonsinglet-quark distribution 
xAg^^, a singlet-quark distribution xAg^, and a gluon distribution xlS.g. As an exam- 
ple, the GS distributions are provided. 
The nonsinglet one is 



QNS<^ =xAu^ + xAd, 



•V 



0.918 * 1.365x°-^^^(l - x)^-^^(l + 11.65X - 4.6v^) 
+ (-0.339) * 3.849x°-^^(l - xf-^^{l + T.Slx - 3.48v^) , 



(5.1) 



the singlet one in the three flavor case is 



QSQ =xAuy + xAdy + QxAS 

=0.918 * 1.365x°-^^2(l - xf-'^^il + 11.65X - 4.6v^) 
+ (-0.339) * 3.849a;°-^^(l - x)^-^^(l + 7.81x - 3.48v^) 
+ 6 * (-0.06) * 18.521x°-^24(^^ _ a,)i4.4(^^ ^ 4 gg^ _ 4 95^) ^ 



(5.2) 



and the gluon distribution is 



GO = 1.71 * 3.099x°-^2^(l - xf-'^^{l + O.Ox + O.Ov^) . 



(5.3) 
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The function QIO calculates two initial flavor distributions. As an example, the GS 
a:Ad+ = xAc? + xAd and xAs~^ — xAs + xAs distributions are provided: 

QIOi ^xAdy + 2xAd 

= - 0.339 * 3.849x°-^^(l - x)^-^%l + 7.81x - 3A8y/x) 

+ 2 * (-0.06) * 18.521x°-^2^(l - xY'^-^l + 4.63x - 4.96v^) , (5.4) 

=2 * (-0.06) * 18.521a;°-^2^(l - x)^^-^(l + 4.63a; - 4.96v^) . (5.5) 
In these equations, the flavor symmetry Au — Ad — As — AS is assumed. 
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6 Numerical analysis 



6.1 Accuracy of evolution results 

There are two parameters which determine the accuracy of our numerical results. They 
are the numbers of steps: and Nf As they become larger, the accuracy becomes 
better. However, there is a restriction because the evolution computation should be 
done within a certain CPU time depending on the machine power. We check how the 
evolution results depend on these parameters and how long it takes for finishing the 
evolution. 

First, we fix the x step at A^2^=1000, then the t step Nt is varied as 20, 50, 200, and 
1000. Numerical results for the polarized singlet quark distribution are shown in Fig. 
[1|. The initial distribution is the GS-A at 0^=^ GeV^ 0. The three-flavor distribution 
is considered. Because there is little information on each antiquark distribution 
flavor symmetric distributions are assumed for the antiquark distributions. The scale 
parameter is taken as A=231 MeV in fitting polarized experimental data. We use 
N f =4 and A=231 MeV for calculating the evolution. The GS-A singlet distribution 
is evolved to the one at (5^=200 GeV^ by our NLO program. The input distribution 
is supplied in the end of the program as the subroutine QSO. The input parameters 
are I0UT=3, IREAD=1, INDIST=1, I0RDER=2, ITYPE=2, IM0RP=2, Q02=4.0, 
Q2=200.0, DLAM=0.231, NF=4, XX=0.0, NX=1000, NT=200, NSTEP=100, and 
NXMIN=— 4. There is almost no difference between the various results in Fig. [1], 
which means that merely fifty steps are enough for getting accurate evolution. This is 
certainly expected from the fact that the scaling violation is a small logarithmic effect. 

Next, the step is varied with fixed A/'t=200 in Fig. |. N^=im, 300, 1000, and 
4000 are taken. Even in the 300 steps, we obtain rather accurate results. The evolution 
results become slightly worse in the gluon case as shown in Figs. |^ and ^. Defining the 
evolution accuracy by | [Aj9(iV^, iV^) - Ap(iV^ = 4000, Nt = 1000)] /Ap{N^ = 4000, A^"* = 
1000) I, we find that the accuracy is better than 1% in the region 10~^ < x < 0.8 with 

= 1000 and A^^ = 200. However, it takes rather a significant amount of running 
CPU time. It is typically seven minutes on the AlphaServer 2100 4/200 in a nonsinglet 
case and sixty minutes in a singlet case. Our program is good enough for a single 
evolution calculation, but it is not efficient enough for repeated use. We will work on 
a much faster program in future. 

We also check our results with other publications. First, the evolution in Ref. 
[0 from (5^=4 GeV^ to Q^=50 GeV^ are compared with our results. The structure 
functions xg^ and xg^ at Q^=50 GeV^ agree well with our evolved ones. Second, the 



evolution in Ref. |]T0[ is compared with ours. The GS-A input distributions are used. 
They are evolved to the ones at Q^=10 GeV^ with A=239 MeV and are also devolved 
to the ones at Q^=2 GeV^. The singlet distribution Ag^ and the gluon distribution 
Ag are compared. We find that they agree well with our evolution results. From 
these studies on our results in comparison with others and from repeated checks on 
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our evolution program, we find that the program BFPl is rehable. We should mention 
that there is another work |TT] on the NLO evolution in addition to Refs. |0 and [ pi| . 



6.2 Calculation of gi from output data 

We do not supply the BFPl program so that the structure function gi is obtained 
directly in the output file. It is because gi depends on the number of flavor and 
because the precise evolution should be described by setting flavor thresholds. We ask 
the reader to set up these points by oneself with our program. Actual evolution results 
for gi are shown in the next subsection, where the following prescription is used for 
calculating gi from the output data file #11. 

First, I0UT=5 or 6 should be chosen so that the evolution of Ag/'(x), Ag^(x), 
Aq^{x), and Ag{x) is calculated by Eqs. ( |2.(j| ) and (f^.Tj). Furthermore, ITYPE=1 
should be chosen for calculating the convolution integrals of these evolved distributions 



with the coefficient functions in Eq. (|2.14|) . The convolution results are written in 



the output file #11, and we denote them as gt^ai^), gi g^ri^)^ fi'is(^)) 9i,gi^)- It 
goes without saying that the convolution is not necessary in the LO case. The LO 
gi is calculated directly by Eq. (p.l3| ). Next, the gi for the proton is, for example, 
calculated by choosing Ad'^{x) and As^{x) for the distributions Aq^{x) and Ag^(x): 

Nf 
i 

1 1 

= - [Agt^ix,Q') ~3{gt/x,Q') + gl,MQ')}] + - Y^e} g,^,{x,Q') , 

(6.1) 

in the three and four flavor cases. Five and six flavor evolution can be calculated in 
the similar way with a slightly modified equation. 



6.3 Comparisons with experimental data 

We calculate evolution of the structure function gi and compare its results with ex- 
perimental data. Using the procedure in the last subsection, we obtain the evolution 
results for gi. The LO and NLO results are shown in Fig. ^. The same GS-A polarized 
parton distributions are used as the input distributions at (5^=4 GeV^ in both LO and 
NLO cases. To be precise, this is not a correct procedure because the GS-A is set up 
for the NLO MS calculation. It should not be used in the LO calculation. Neverthe- 
less, the GS-A distributions are also employed in the LO for finding NLO effects on 
gi. Therefore, the differences between the solid and dashed curves at (5^=4 GeV^ are 
purely due to the coefficient functions. Because the coefficient-function contributions 
are rather large, it is difficult to find the NLO evolution effects from Fig. |^. In order to 
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illustrate the NLO effects from Q^=4 GeV^ to (^^=200 GeV^, we show the evolution 
of the distribution in Fig. |^. The same distributions are taken at Q^=4: 

GeV^ in the LO and NLO cases. The figure indicates that the NLO effects are posi- 
tive in the small and large x regions and they are negative in the intermediate region 
X ^ 0.2. 

Our evolution results are compared with spin-asymmetry data in Figs. ||, and 
|. The SLAC-E130 |T2|, SLAC-E143 0, EMC 0, and SMC data are shown 



in the figures. We comment how the spin asymmetry is calculated theoretically. The 
asymmetry is given by the structure functions gi and Fi as 

The function R is given by i? = {F2 — 2xFi)/{2xFi), and it could be taken from 



Ref. |jT6|. In fitting experimental data for obtaining the optimum unpolarized parton 
distributions, the F2 structure function is used instead of Fi. Therefore, it is better to 
calculate the asymmetry with F2 if we would like to compare with the experimental 
data. The MRS-G unpolarized parton distributions [O] are used for calculating F2 



and the function R in Ref. [T^ is used in Eq. (|6.2|). 

In Fig. 1^, our evolution curves at a;=0.035 are shown with the asymmetry Ai data. 
The dashed and solid curves indicate the LO and NLO evolution results respectively. In 
the large region, both results are almost the same; however, they differ significantly 
at small in particular in the region <2 GeV^. The difference is not so large at 
slightly larger x (=0.08) as shown in Fig. |^. However, the LO values become larger than 
those of the NLO evolution. In the medium x region (a;=0.25), the difference becomes 
larger again at small as shown in in Fig. ^ From these figures, we find that the 
asymmetry has dependence although it is not large. People used to assume that the 
asymmetry is independent of by neglecting the evolution difference between gi 
and Fi in analyzing the experimental data. We find clearly that it is not the case. For 
a precise analysis, the dependence in the asymmetry has to be taken into account. 
Currently, parametrization studies are in progress by using our program. We expect 
to obtain NLO fitting results in the near future. 
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7 Summary 



We have investigated numerical solution of the spin-dependent DGLAP equations with 
or without the NLO corrections. The solution method is so called brute-force method. 
A FORTRAN program is provided for evolving longitudinally polarized parton dis- 
tributions including nonsinglet, singlet, each flavor, and gluon distributions. Further- 
more, the evolution results could be written in a structure-function form, which is the 
convolution of the evolved distributions with the coefficient functions. Therefore, the 
structure function gi could be also calculated from our output data. We checked that 
typical accuracy is better than 1% with reasonable numbers of steps: A^a;=1000 and 
A^j=200. Comparisons with experimental data indicate significant dependence and 
NLO effects in the region Q"^ k, 1 GeV^. Therefore, the independence assumption, 
which was used to be employed in analyzing Ai experimental data, is not valid. We 
have to include the evolution differences between Fi and gi. Our evolution code is 
very useful for theoretical and experimental researchers in high-energy spin physics. 
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Appendix A. Splitting functions 

First, the constants Pq, /?i Cg, Cp, and Tr are given by the number of color {Nc=3) 
and the number of flavor (Nf) as 



(A.l) 



G 



2Nr 



R 



Splitting function in the leading order(LO) are 



iTiJ 
AP^f (x) 



(l-x)^ 
5ij AP^l{x) , 

2Tr{2x-1) 

Cf {2-x) , 
1 



1 - X + - 6{1 - x) 



2C, 



a 



(l-x). 



-2x + l 



+ ^Sil-x) 



where the + function is defined by 











1 — X 



(A.2) 



(A.3) 



(A.4) 



It should be noted that the above integration is defined in the region < x < 1. 
The NLO splitting functions for Aqf{x) and Ag{x) are given by 



AP 



(1) 



x] 



<ii y 

AP^^lix) 

9q7 



CFTnF;^ix) + CGTnFl{x) , 

Cf Tn Nj FUx) + Cl p-'ix) + Cf Cg F'ix) 



(A.5) 



AP^l\x) = -Cg Tn Nf F^^ix) - Cf Tr Nj F^^x) + C^ F^^x) . 
The splitting functions in the singlet equations ( p^ ) are expressed as 

APqq ^Aqs = J2 ^^44 ® ' ^^99 ® Ag = ^P.+g ® ^9 

APgg ® Ag, = ^ AP^^+ ® Aqj . 



(A.6) 
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The splitting function APj± is given 



A p(l) _ p(l) 



(A.7) 



where j^fg are the unpolarized NLO nonsinglet sphtting functions. The ± in these 
equations indicates "Ag ± Aq type" distribution ^^aj(Agj ± Agj). The function 



^^g± TVS is giv6^ by 



AP^(i)^^(x) <! Pp{x) T Pa{x) + 5(1 - x) 



^-^^' + C(3)- 8^(00) 



+ -CpC^<{ Pg{x) ± P^(x) + 5(1 - x) 



^ + y^'-C(3) + 8S(oo) 



+ Ci.TRiV^<| PnAx) - 5(1 - x) ( i + ^tt^ 



(A.8) 



where Pf{x), Pg{x), Pnf{x), and Pa (a;) are given in Ref. 
Pf(x) 
Pg{x) 



2 Inx ln(l — a;) — I \- 2x ] Inx (1 + a;) In^ a; —5(1 — x) , 

V 1 — X / 2 



1 — X 
1 + 



(1-,t)_ 



11 

y 



67 1 



In X -\ lnx + -: -TT 



Pnp{x) 
Pa{x) 



1+x^ 



In X 

3 



(l-x)_ 

l + a;2 fm+^)dz^ l-z 
— m 



9 3 
2(1 -x) 



40 



+ 2(l + x) lnx + — (1-x) 



+ 2(1 + x) Inx + 4(1 -x) . 



1 + X Jx/{l+x) Z 

The functions Fqq, Fgg, Fgg, and Fgg are defined by 

Fqq{x) = (1 — x) — (1 — 3x) Inx — (1 + x) In^ x 



(A.9) 



(A.IO) 



22 + 27x- 91nx + 8(l -x) ln(l - x) 



+ 5pgg{x) 



21n^(l -x) -41n(l -x)lnx + ln^x tt^ 



2(12 - llx) - 8(1 - x) ln(l - x) + 2(1 + 8x) Inx 



IT 



ln^(l -x) 

^ ' 6 



Spqg{x) - [2S'2(x) -31n^x] 6pqg{-x) 



(A.ll) 
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Fgq{x) = - ^ - ^(4 - 2;) Inx - Spgg{-x) ln(l - x) + 



-4 - ln^(l - x) + -Wx 



F^g{x) = (4 - 13a;) Inx + -(10 + x) ln(l - a;) + -(41 + 35a;) 



+ -[-2^2(2;) + Sin' x]5pgg{-x) 



+ 



ln'(l — x) — 21n(l — x)lnx — 



9 



Spgg{x) , 



(A.12) 



4 20 4 

^i(2;) =4(1 - a;) + -(1 + x) Ina; + y%,(a:) + -S{1 - x) , 

F^g{x) =10(1 - a;) +2(5-a;)lna; + 2(l + a;)ln'a; + (5(l -x) , 

1 19 
Fgg{x) = -(29 - 67a;) Ina; - y (1 - a;) + 4(1 + a;) In' a; - 2S2{x)6pgg{-x) 



3 
+ 



— - 4 ln(l - a;) In a; + W a; - — 
9 o 



5pgg(a;) + 



3C(3) + 



where, Spgg , Spgg , Spgg Rie defined 



Spqg{x) = 2X — 1 

^Pgq{x) = 2 - a; , 
1 

%g(x) 



5(1 - X) , 

(A.13) 
(A. 14) 



2a; + 1 . 



The C function is defined by C(/i;) = X^^i C(3) is given by the numerical value 

(C(3)=l. 2020569...). The S2 function is expressed in terms of the Spence function S{x) 

dz^l-z 

Jx/{l+x) ^ ^ 

X 



Soix) 



S 



1 + x 

where S{x) is defined by 



S 



1 + x 



In' 



1 + a; 



In' 



X 



1 + x 



S{x) 



^ , Inz 
dz 



z 



(A.15) 



(A.16) 



It should be noted that another convention is sometimes used, namely —S{x) may be 
called the Spence function. It is useful to use a series expansion form for numerical 
calculations 

(1-a;)*^ 



fe=i 



A;2 



(A.17) 
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Appendix B. Coefficient functions 

The coefficient functions ACg and ACg are 



In 



ACg{x,as 



2^ 



ABJx 



(B.18) 
(B.19) 



where the functions ABg and ABg are 



ABgix) 



{l + x') 



2. I ln(l - x) 



3 1 



+ 2 + X 



1 + a;^ 

1 ~ X j, 2(1 — x)+ 1 — X 



5 + 



In a; 



(B.20) 



(2x - 1) (in 



1 — X 



\ X 



1 +2(1 -x) 



(B.21) 
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TEST RUN OUTPUT 



IOUT= 3 IREAD= 1 INDIST= 1 IORDER= 2 ITYPE= 2 IMORP= 2 

Q02= 4.0000 Q2= 200.000 DLAM= 0.2310 NF= 4 

XX= 0.0000000 NX=1000 NT= 200 NSTEP= 50 NXMIN= -4 



0.000100 
0.000120 
0.000145 
0.000174 
0.000209 
0.000251 
0.000302 
0.000363 
0.000437 
0.000525 
0.000631 
0.000759 
0.000912 
0.001096 
0.001318 
0.001585 
0.001905 
0.002291 
0.002754 
0.003311 
0.003981 
0.004786 
0.005754 
0.006918 
0.008318 
0.010000 
0.012023 
0.014454 
0.017378 
0.020893 
0.025119 
0.030200 
0.036308 
0.043652 
0.052481 
0.063096 



-0.019166 

-0.020672 

-0.022276 

-0.023982 

-0.025790 

-0.027703 

-0.029718 

-0.031833 

-0.034043 

-0.036338 

-0.038707 

-0.041130 

-0.043584 

-0.046036 

-0.048444 

-0.050754 

-0.052898 

-0.054790 

-0.056324 

-0.057372 

-0.057777 

-0.057353 

-0.055883 

-0.053115 

-0.048770 

-0.042545 

-0.034126 

-0.023214 

-0.009552 

0.007027 

0.026539 

0.048784 

0.073276 

0.099181 

0.125309 

0.150159 



0.097380 
0.105178 
0.113534 
0.122478 
0.132045 
0.142266 
0.153174 
0.164803 
0.177185 
0.190351 
0.204330 
0.219150 
0.234833 
0.251398 
0.268856 
0.287209 
0.306451 
0.326559 
0.347493 
0.369192 
0.391566 
0.414496 
0.437819 
0.461325 
0.484750 
0.507761 
0.529949 
0.550821 
0.569787 
0.586160 
0.599147 
0.607861 
0.611331 
0.608536 
0.598449 
0.580109 
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0.075858 
0.091201 
0.109648 
0.131826 
0.158489 
0.190546 
0.229087 
0.275423 
0.331131 
0.398107 
0.478630 
0.575440 
0.691831 
0.831764 
1.000000 



0.172073 
0.189500 
0.201326 
0.207180 
0.207504 
0.203229 
0.194990 
0.182213 
0.162801 
0.134220 
0.096195 
0.053887 
0.018721 
0.001943 
0.000000 



0.552726 
0.515809 
0.469336 
0.413952 
0.351177 
0.283585 
0.214868 
0.149663 
0.092992 
0.049196 
0.020506 
0.005833 
0.000825 
0.000022 
0.000000 



Figure 1: Nt dependence of singlet evolution results is shown. The initial distribution 
is the GS-A at Q2=4 GeV^. It is evolved to the one at Q2=200 GeV^ by the next-to- 
leading-order DGLAP evolution equations. A^j;=1000 is fixed and Nf is varied (A^t=20, 
50, 200, and 1000). There are dotted, dashed, dot-dashed, and solid curves at Q'^—200 
GeV^ for A^t=20, 50, 200, and 1000 respectively. 
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Figure 2: N^. dependence of singlet evolution results is shown. The initial distribution 
is the same in Fig. |I]. Nt=200 is fixed and is varied (iV^=100, 300, 1000, and 4000). 
There are dotted, dashed, dot-dashed, and solid curves at Q^=200 GeV^ for iV3;=100, 
300, 1000, and 4000 respectively. 
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Figure 5: evolution results for the proton structure function g\ are shown. The 
solid curves are xg\ at Q'^—^ and 200 GeV^ in the NLO case, and the dashed ones are 
in the LO. 
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Figure 6: evolution results for the distribution are shown. The same 

distributions are assumed at Q^=4 GeV^ and the evolved distributions are shown by 
the solid curve in the NLO case and by the dashed one in the LO at 200 GeV^. 
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Figure 7: dependence of the spin asymmetry for the proton is calculated in the 
LO and NLO cases. Calculated results at x = 0.035 are compared with SLAC-143 [|I^, 
and SMC [jl5| experimental data. The solid curve is in the NLO case and 
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the dashed one is in the LO. 
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Figure 8: dependence of the spin asymmetry at x = 0.08. Notations are same 
with those in Fig. 0. 
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Figure 9: dependence of the spin asymmetry at x = 0.25. Notations are same 
with those in Fig. % The SLAC-E130 data [|T2[ is added. 
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